% =========================================================================
% modeFinderInitialize
% Generates starting parameter values
% Initializes some of the loop_ matrices 
% A. Justiniano Feb 21 2014 
% 
% Function needs to be cleaned a lot, once incorporate priors and parameters 
% in structures. 
% Still numerous variables created that are not used elsewhere. 
% =========================================================================
%% Uses modePosterior.m for the optimization 
posStru.param.cal=posStru.param.cal; 
posStru.steady.est=ssposest; 

parStru.param=zeros(numpar, 1);
parStru.param(posStru.param.cal)=parStru.cal; 

dataStru.data=Y; 
dataStru.trainVec=trainvec; 

filterStru.tauVec=addsol.tauVec; 
filterStru.funcLikel=addsol.funcLikel; 
filterStru.funcKfilter=addsol.funcKfilter; 
filterStru.funcPost=@modePosterior; 
filterStru.funcMin=@modePosteriorMin; 

offset=1e-7; 


%% posStru.param.estimated posStru.param.est 
%% pos
%%            
display('________________________________'); 
display('Begin Generating Starting Values'); 

npoints_full = 0; 
if priorStru.randomGridStart==0 && priorStru.randomPriorStart==0;
    % use loopStru.parStartMat
    disp('Using Starting values');    
    loopStru.parStartMat=parstrvals; 
    priorStru.numStartVals=size(loopStru.parStartMat,2); 
else   
    loopStru.parStartMat= zeros(numpar,priorStru.numStartVals); 
    icount = 0; % iteration counter
    fcount = 0; % fill count;
    disp('_________________________________________________');
    disp('Starting values, uniform from grid');
    pause(0.25);
    
    disp('_________________________________________________');
    disp('Starting values, uniform from prior');
    pause(0.25);
    
    while npoints_full == 0;
        admissible = 0;
        icount = icount + 1;
        if icount > 100000;
            error('100000 draws cannot generate admissable value')
        end
        parStru.param = zeros(numpar, 1);
        switch priorStru.randomGridStart==1
            case 1 % random from uniform
                for ii = 1:numpar;
                    if  pargrid(ii, 2) - pargrid(ii, 1) > 2*offset
                        parStru.param(ii) = unifrnd(pargrid(ii, 1)+offset, pargrid(ii, 2)-offset);
                    else
                        parStru.param(ii) = unifrnd(pargrid(ii, 1), pargrid(ii, 2));
                    end
                end
                parStru.param(posStru.param.cal) = parStru.cal;
            case 0
                % random from prior
                tempdist = char(prpar.dist(posStru.param.est));
                parstrval_est_temp = genmultrand(tempdist,prpar.alphap(posStru.param.est),...
                    prpar.betap(posStru.param.est),1,prpar.lbnd(posStru.param.est),...
                    prpar.ubnd(posStru.param.est));
                parStru.param(posStru.param.est) = parstrval_est_temp;
                parStru.param(posStru.param.cal) = parStru.cal;
        end
        
        tone=parStru.param(posStru.param.est)-prpar.lbnd(posStru.param.est);
        two=prpar.ubnd(posStru.param.est) - parStru.param(posStru.param.est);
        if all(tone >= 0) == 1 && all(two >= 0) == 1;
            % Evaluate likelihood and posterior            
            
%             [lpostd,likel,lpriord,lprior_ss,flag_ok]=feval(@lmjPosterior,...
%                 parStru.param(posStru.param.est),parStru.param,...
%                 posStru.param.est,model.handle, Y, trainvec, prpar, prss, solveopt, addsol, ssposest); 
%            
            
        [lpostd,likel,lpriord,lprior_ss,flag_ok]=feval(filterStru.funcPost,...
            parStru.param(posStru.param.est),parStru,...
        posStru,model,dataStru,prpar,prss,filterStru,flags);     

            disp(' ');
            disp(sprintf('LogLikel=%10.2f',likel));
            disp(sprintf('lpriord=%10.2f',lpriord));
            disp(sprintf('lprior_ss=%10.2f',lprior_ss));
            disp(sprintf('LogPosterior=%5.5e',lpostd));
            
            if lpostd > -1e15 && flag_ok==1
                admissible = 1; 
            end 
        end 
            if admissible == 1;
                fcount = fcount + 1; % Keep or not
                loopStru.parStartMat(:, fcount) = parStru.param; 
            end
            npoints_full = (fcount == priorStru.numStartVals);  
    end 
end
%%

clear fcount npoints_full admissible lpostd two tone loop_logPost offset;  
clear temp*; 


%% loop_ functions initilized here because parfor does not allow double indexing i.e. 
%  loop_logPrior(ii,2) and loop_logPrior(ii,1)
%  Also parfor does not like it in structures 
% =========================================================================
disp('Checking initial densities....');
loop_logPrior(:,2)=-inf*ones(priorStru.numStartVals,1);
loop_logLikel(:,2)=-inf*ones(priorStru.numStartVals,1);
loop_logPost(:,2)=-inf*ones(priorStru.numStartVals,1);
for ii = 1:priorStru.numStartVals;
    parStru.param=loopStru.parStartMat(:, ii);
    [loop_logPost(ii,1),loop_logLikel(ii,1),loop_logPrior(ii, 1)]=...
        feval(filterStru.funcPost,loopStru.parStartMat(posStru.param.est,ii),parStru,...
        posStru,model,dataStru,prpar,prss,filterStru,flags);     
    if loop_logPost(ii,1) < -1e15
        error('Likelihood too small! Parameter value is not admissible')
    end
    dispaj('Starting densities for candidate ',ii);
    disp([sprintf('LogPosterior=%5.5e',loop_logPost(ii,1)),sprintf('LogLikelihood=%5.5e',loop_logLikel(ii,1))]);   
    if abs( loop_logPrior(ii,1) -( loop_logPost(ii,1) - loop_logLikel(ii,1) ) ) > 1e-3;
        error('Log Posterior- Log Likelihood does not equal Log Posterior') 
    end 
end
disp('Obtained PARSTRVALS');
savecell([parnames num2cprec(loopStru.parStartMat,7)],[],location.savePath,'starting values');

